module mod_chem_transport

   use iso_fortran_env, only: dp => real64
   implicit none
   private
   public :: type_chem_transport, type_chem_gas_transport

   type :: type_chem_transport
      character(len=:), allocatable :: name
      character(len=:), allocatable :: geometry
      real(dp), allocatable :: viscosity_coefficients(:)
      real(dp), allocatable :: thermal_conductivity_coefficients(:)
   contains
      procedure :: get_viscosity => GetViscosity
      procedure :: get_thermal_conductivity => GetThermalConductivity

   end type type_chem_transport

   type :: type_chem_gas_transport
      character(len=:), allocatable :: name

   contains
      procedure :: get_viscosity => IdealGasGetViscosity
      procedure :: get_thermal_conductivity => IdealGasGetThermalConductivity
      procedure :: get_phi => IdealGasGetPhi
   end type type_chem_gas_transport

contains

   pure function GetViscosity(this, T) result(viscosity)
      class(type_chem_transport), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: viscosity
      integer :: l, n

      n = size(this%viscosity_coefficients)
      viscosity = exp(sum([(this%viscosity_coefficients(l)*log(T)**(l - 1), l=1, n)]))

   end function GetViscosity

   pure function GetThermalConductivity(this, T) result(thermal_conductivity)
      class(type_chem_transport), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: thermal_conductivity
      integer :: l, n

      n = size(this%thermal_conductivity_coefficients)
      thermal_conductivity = exp(sum([(this%thermal_conductivity_coefficients(l)*log(T)**(l - 1), l=1, n)]))

   end function GetThermalConductivity

   function IdealGasGetViscosity(this, X, molar_weights, partial_viscosity) result(viscosity)
      class(type_chem_gas_transport), intent(in) :: this
      real(dp), intent(in) :: X(:)
      real(dp), intent(in) :: molar_weights(:)
      real(dp), intent(in) :: partial_viscosity(:)
      real(dp) :: viscosity

      real(dp), allocatable :: phi(:, :)
      integer :: i

      phi = this%get_phi(molar_weights, partial_viscosity)
      viscosity = 0.0_dp

      do i = 1, size(X)
         viscosity = viscosity + X(i)*partial_viscosity(i)/sum(X*phi(i, :))
      end do

   end function IdealGasGetViscosity

   function IdealGasGetThermalConductivity(this, X, molar_weights, partial_thermal_conductivity) result(thermal_conductivity)
      class(type_chem_gas_transport), intent(in) :: this
      real(dp), intent(in) :: X(:)
      real(dp), intent(in) :: molar_weights(:)
      real(dp), intent(in) :: partial_thermal_conductivity(:)
      real(dp) :: thermal_conductivity

      ! real(dp), allocatable :: phi(:, :)
      ! integer :: i

      ! phi = this%get_phi(molar_weights, partial_thermal_conductivity)
      ! thermal_conductivity = 0.0_dp
      ! do i = 1, size(X)
      !    thermal_conductivity = thermal_conductivity + X(i)*partial_thermal_conductivity(i)/sum(X*phi(i, :))
      ! end do
      thermal_conductivity = 0.5_dp*(sum(X*partial_thermal_conductivity) + 1.0_dp/sum(X/partial_thermal_conductivity))

   end function IdealGasGetThermalConductivity

   function IdealGasGetPhi(this, M, mu) result(phi)
      class(type_chem_gas_transport), intent(in) :: this
      real(dp), intent(in) :: M(:)
      real(dp), intent(in) :: mu(:)
      real(dp), allocatable :: phi(:, :)

      integer :: i, j, n
      real(dp), parameter :: rROOT8 = 1.0_dp/sqrt(8.0_dp)

      n = size(M)

      allocate (phi(n, n))

      do j = 1, n
         do i = 1, n
            phi(i, j) = rROOT8*(1.0_dp + M(i)/M(j))**(-0.5_dp)*(1.0_dp + sqrt(mu(i)/mu(j))*(M(j)/M(i))**0.25_dp)**(2.0_dp)
         end do
      end do

   end function IdealGasGetPhi

end module mod_chem_transport
